ADA104578 


UNCLASSIFIED _ 

SECURITY  CLASSIFICATION  OR  THIS  PACE  (When  Oaf  a  Entered) 


1  REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

mm+rmi  number 

..  .aCL-.TR- 19  4  4 _ A 

2.  GOVT  ACCESSION  NO. 

Ab~A 

3.  RECIPIENT'S  CATALOG  NUMBER 

L£Z£ _ 

U.  nrtf  fw  «mm wta) 

|/  The  Deconvolution  of  Aerosol  Back- 
scattered  Optical  Pulses  to  Obtain 
System- Independent  Aerosol  Signatures  . 


w 

Technical  Report ■ 


rYPE  OF  REPORT  ft  PERIOD  COVEREO 


6.  PERFORMING  ORG.  REPORT  NUMBER 


7.  AUThOR/AJ 

8.  CONTRACT  OR  GRANT  NUMBER/*.) 

Dennis  McGuire 

Michael,'  Conner 

S.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

10.  PROGRAM  ELEMENT.  PROJECT.  TASK 

Harry  Diamond  Laboratories 

AREA  ft  WORK  UNIT  NUMBERS 

2800  Powder  Mill  Road 

Adelphi,  MD  20783 

Program  Ele. :  MAIF 

11.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

12.  REPORT  DATE 

Commander  ,  ■  ) 

U .  S .  Army  ARRADCOM  '  * 

June  1981 

.  U.  NUMBER  OF- PAGES  ] 

Dover,  NJ  07801 

66 

14.  MONITORING  AGENCY  NAME  ft  ADORESSf/f  different  from  Controlling  OWca ) 

IS.  SECURITY  CLASS,  (of  thle  report ) 

16-  DISTRIBUTION  STATEMENT  (of  thl«  Report) 


UNCLASSIFIED _ 

15*.  OECL  ASSl  FI  CATION/  DOWNGRADING 
SCHEDULE 


Approved  for  public  release;  distribution  unlimited. 


17.  DISTRIBUTION  STATEMENT  fof  Ida  adafracf  anfar.rf  In  Stock  20,  It  dttteeent  from  Report) 


it.  supplementary  notes 

HDL  Project:  A18032 
DRCMS  Code:  36KA5000204 


IS.  KEY  WOROS  (Continue  on  favaraa  aid*  It  nacaaaary  and  Idantlty  by  block  numbar) 

Aerosol  backscatter  Optical  fuzing 

Deconvolution  Active  optical  sensors 


20.  ABSTRACT  (CantBuu*  aaa  rararma  at*a  If  nacaaaary  mtd  I  dan/ tty  by  block  numbar) 

Means  are  discussed  for  extracting  system-independent  aerosol 
signatures  from  aerosol  backscatter  measurements  obtained  with  a 
specific  pencil  beam  active  optical  detection  system.  Such 
signatures  are  required  before  the  backscatter  data  can  be  applied 
to  various  proposed  optical  fuze  designs  for  determining  their 
aerosol  vulnerability  and  to  the  investigation  of  aerosol 
discrimination  schemes.  The  measurement  system,  which  has  been  — 


00  ,12^1  1473  eo.no«  or  *  novo  is  obsolete  UNCLASSIFIED 

SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Dole  Entered) 


1 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAOE(TFh*n  DMm  Snl.mtf) 


ABSTRACT  (Cont'd) 

used  in  numerous  experiments  to  probe  such  aerosols  as 
weather  clouds  and  military  smokes,  is  a  short  pulse  GaAs 
laser  probe  (pulse  width  sio  ns)  whose  range  sensitivity 
extends  from  near  the  system  to  beyond  10  m. 

The  resolution  needed  for  realistic  evaluations  of  aerosol 
density  gradient  effects  on  optical  fuzing  systems  is  greater 
than  that  inherent  in  the  spatial  width  of  the  probing  pulses, 
which  is  approximately  2  m.  To  achieve  such  superresolution 
requires  analytical  and  numerical  techniques  to  deconvolve 
the  transmitted  pulses  from  the  signatures.  Deconvolution 
methods  are  investigated  for  measurement  noise  conditions 
typical  of  the  backscatter  data .1  A  computationally  fast 
numerical  deconvolution  algorithm  is  devised  together  with  a 
comprehensive  supporting  analysis.  Both  indicate  that  severe 
signal-to-noise  ratio  constraints  apply  to  the  achievement  of 
meaningful  superresolution. 

While  the  signal-to-noise  ratios  typical  of  recent 
measurements  are  likely  to  satisfy  the  severe  constraints 
discovered,  many  of  the  earlier  data  are  too  noisy  and  thus 
require  other  signature  determination  methods^  An  alternative 
approach  is  discussed  and  given  a  preliminary^bvaluation. 

The  approach,  which  seems  promising,  involves  the  use  of  a 
priori  knowledge  of  the  nature  of  the  signature  coupled  with 
parameter  estimation  techniques. 


eoesslon  For 

NTIS  GRA&I 
DTIC  TAB  □ 

Unannounced  □ 

Just  if  ication__ _ 

By _ _ 

Dist  ri but i on/ 
Availability  Codas 
Avail  and/or 
1st  Special 


ELECTE 

SEP  2  0  1501 


RE:  Classified  References,  Distribution 
Unlimited 

No  change  in  distribution  statement  per 
Mr.  Dennis  McGuire,  HDL/DELHD-RT-CB 


CONTENTS 


1.  INTRODUCTION . 7 

2.  PRELIMINARY  DISCUSSION  OF  SIGNATURE  EXTRACTION  PROBLEM . 11 

2.1  Linear  System  Description... . 11 

2.2  Typical  Transmitter  Pulses,  Range  Law,  and  Signatures . 12 

3.  FORMAL  AND  NUMERICAL  SOLUTIONS  OF  DECONVOLUTION  PROBLEM . 20 

3.1  Formal  Solution  and  Effect  of  Noise . 20 

3.2  Numerical  Deconvolution . 23 

3.2.1  Discrete  Formulation . 23 

3.2.2  Z-Transform  Deconvolution . 25 

4.  EFFECT  OF  NOISE  IN  SIGNATURE  EXTRACTION . 31 

4.1  Numerical  Effect. . 31 

4.2  Analytical  Estimation  of  Noise  Effects . 33 

5.  ALTERNATIVE  APPROACH  THROUGH  PARAMETER  ESTIMATION . 44 

5.1  Estimation  of  o  and  w  for  Uniform  Aerosol . 46 

5.2  Estimation  of  o  and  M  from  Measured  Return  Signals . 51 

6.  SUMMARY  AND  DISCUSSION . 60 

LITERATURE  CITED . 63 

DISTRIBUTION . 65 


3 


FIGURES 


1  Shape  of  typical  optical  pulse  from  GaAs  laser  pulser . 13 

2  Shape  of  typical  optical  pulse  from  improved  GaAs  laser 

pulser. . . ....14 

3  Measured  range-response  curve,  R(x),  of  pencil  beam  GaAs  laser 

probe  for  aerosol  backscatter  measurements . 15 

4  Illustration  of  idealized  aerosol  density  profile  along  pencil 

beam  influence  pattern  of  active  optical  detection  system . 17 

5  Aerosol  signatures  calculated  from  defining  equation  (1)  for 

aerosol  density  profile . 18 

6  Aerosol  signature  for  semi-infinite  uniform  aerosol 

di  str  ibution . 19 

7  Computed  aerosol  signature  according  to  equations  (6) 

and  (7) . 27 

8  Product  of  aerosol  signature  of  figure  7  and  range-response 

function  computed  from  equation  (29) . . . 28 

9  Computed  aerosol  return  signal  according  to  equation  (23) . 29 

10  Product  of  aerosol  signature  and  range  response  function . 30 

1 1  Result  of  deconvolving  return  signal  of  figure  9.. . 32 

12  Comparison  of  shape  of  model  transmitter  pulse  with  that  of 

measured  transmitter  pulse  of  figure  2 . 34 

13  Sketch  of  graph  of  I(x)  =  (^x)^(l  -  x^)^/sin2  irx . ...36 

x 

/Q 

I(x)  dx  versus  xc  and  versus 

fc  =  xc/2T . 7 . 38 

15  Signal-to-noise  ratio  penalty  factor  for  model  transmitter 

pulse  of  equation  (32) . . . . . 42 

16  Signal-to-noise  ratio  penalty  factor  for  measured  transmitter 

pulse  of  figure  2............. . 43 


4 


FIGURES  (Cont'd) 


Constant-cost  contours  of  J(o,u) . 

Typical  trajectories  for  steepest  descent  algorithm... 

Typical  trajectories  for  intuitive  algorithm . 

Sampled  values  of  measured  cumulus  cloud  return  signal 


1 


INTRODUCTION 


A  measurement  program  that  collects  laser  backscatter  data  from 
various  aerosols  and  simultaneously  makes  measurements  to  characterize 
the  aerosol  environment  has  been  underway  at  the  Harry  Diamond  Labora¬ 
tories  ( HDL)  for  several  years.  Measurements  have  been  made  on  water 
clouds  during  helicopter  flight  tests  and  on  other  aerosols  such  as 
smoke,  dust,  and  fog.  The  purpose  of  these  measurements  has  been 
largely  to  furnish  the  data  needed  to  determine  aerosol  backscatter 
effects  on  active  optical  fuzing  (AOF)  systems  using  GaAs  injection 
laser  transmitters,  since  such  effects  could  cause  false  target 
problems.  This  report  discusses  the  results  of  several  related  investi¬ 
gations  whose  purpose  was  to  provide  the  needed  interpretations  of  the 
backscatter  data  and  the  methods  by  which  the  data  can  be  used  to  deter¬ 
mine  aerosol  backscatter  effects  in  various  proposed  AOF  systems.  The 
measurement  program  itself  is  discussed  fully  by  McGuire,  Smalley,  and 
Sztankay . 

The  methods  used  for  backscatter  signal  acquisition  and  subsequent 
digitization  of  the  data  are  discussed  in  detail  by  Vanderwall  and 
Conner.3 4  The  data,  a  train  of  received  signal  pulses,  are  recorded 
during  the  measurements  on  video  tape  in  standard  television  (TV)  format 
so  that  each  TV  frame  contains  the  data  for  an  entire  backscattered 
pulse.  A  specially  designed  electronic  frame-code  generator  applies  a 
sequential  code  number  to  each  TV  frame  during  the  data  recording.  By 
using  another  specially  designed  device,  the  coded  video  tape  data  are 
automatically  digitized  in  our  laboratory  and  thereby  put  in  a  form 
suitable  for  further  processing  and  analysis  by  computer. 

The  immediate  question  for  data  usage  concerns  how  backscatter  data 
obtained  with  a  particular  active  optical  detection  system — the  measure¬ 
ment  system  employs  a  GaAs  short-pulse  laser  probe — can  be  translated  to 
apply  to  another  system  that  might  differ  from  the  first  in  such 
respects  as  range  coverage  and  transmitter  pulse  characteristics.  An 
answer  can  be  formulated  in  terms  of  an  aerosol  signature  concept, 
provided  that  multiple  scattering  contributions  to  the  aerosol  return 
signals  can  be  ignored. 

3D.  McGuire,  H.  M .  Smalley,  and  Z.  G.  Sztankay ,  Measurements  of 
Backscatter  Effects  in  Clouds  at  0.9  vim  (U ) ,  Proc.  JTCG/MD/WPFF  Tri- 
Service  Optical  Fuze  Technology  Symposium,  Naval  Weapons  Center  NWC  TP 
5871,  Part  I  (October  1976),  45-74.  (CONFIDENTIAL) 

2Z.  G.  Sztankay  and  D.  McGuire,  Backscatter  in  Clouds  at  0.9  vim  and 
Its  Effects  on  Optical  Fuzing  Systems  (U),  Proc.  Seventh  DoD  Conference 
on  Laser  Technology  (November  1977).  (SECRET) 

3Z.  G.  Sztankay,  Measurement  of  the  Localized  Optical  Charac¬ 
teristics  of  Natural  Aerosols ,  Smoke,  and  Dust,  Proc.  Smoke/Obscurants 
Symposium  II  (25-26  April  1979). 

4 Jonathan  Vanderwall  and  Michael  Conner,  A  Novel  Scan-Converting 
Oscillographic  Technique  for  In-Situ  Signal  Acquisition  with  Subsequent 
Automatic  Digitization,  Harry  Diamond  Laboratories  HDL-TR-1956  (1981). 

7 


FHLCLUUG  FAuK  bb*J*-HOI 


To  develop  the  aerosol  signature  concept,  consider  first  the  GaAs 
laser  probe  used  to  collect  the  backscatter  data.  It  is  a  pencil  beam 
active  optical  detection  system  employing  a  short-pulse  GaAs  injection 
laser  transmitter  and  a  transceiver  influence  pattern  sensitive  from 
ranges  near  the  system  out  to  ranges  in  excess  of  10  m.  The  transmitter 
pulse  widths,  measured  by  the  full  width  at  half  maximum  (FWHM),  are 
currently  5  ns,  but  many  of  the  data  were  obtained  with  wider  pulses  (7 
and  11  ns  FWHM)  using  earlier  versions  of  the  probe.  The  shape  of  the 
detected  aerosol  return  pulses  depends  on  several  factors  besides  the 
distribution  of  aerosol  in  the  influence  pattern.  These  factors  include 
the  shape  of  the  transmitter  pulses  and  the  shape  of  the  range  sensi¬ 
tivity  curve  of  the  system. 


The  foregoing  dependencies  can  be  expressed  in  a  useful  mathematical 
form  with  the  aid  of  the  following  definitions.  Let  P(t)  and  V(t)  be 
the  instantaneous  transmitted  power  and  return  signal,  respectively,  and 
let  x  denote  range  from  the  transceiver  measured  along  its  pencil  beam 
influence  pattern.  Define  the  function  of  range  C(x)  by 


C(x)  =  y(x)  exp 


o(  s) 


(1) 


where  u(x)  and  o(x)  are,  respectively,  the  volume  backscatter  and 
extinction  coefficients  of  the  aerosol  at  range  x.  Finally,  let  R(x) 
denote  the  range  sensitivity  function  of  the  system.  (The  function  R  is 
defined  as  the  peak  receiver  response  to  a  flat  reflecting  target, 
oriented  at  right  angles  to  and  filling  the  influence  pattern,  as  a 
function  of  range  from  the  transceiver  to  the  target.  The  scale  of 
values  of  R  depends  on  the  units  of  receiver  response  and  the  target 
reflectivity;  ordinarily,  some  convenient  normalization  of  the  scale  is 
used.)  Then  V( t)  can  be  expressed  as 


V(t)  *  K  /“  P(t-x )C(ct/2)R(ct/2)  dt  , 


(2) 


provided  that  the  receivers  respond  linearly  to  the  received  optical 
power;  the  factor  K  is  a  constant  depending  on  the  normalization  chosen 
for  R(x)  and  the  units  of  V(t);  c  is  the  speed  of  light.  This  result, 
originally  derived  by  Burroughs,^  has  considerable  generality;  it 
applies  not  only  to  the  GaAs  laser  probe,  but  to  virtually  all  pencil 


5tf.  H .  Burroughs,  Computation  of  Cloud  Backscatter  Power  as  a 
Function  of  Time  for  an  Active  Optical  Radar  (U ) ,  Naval  Weapons  Center 
NWC  TP  5090  (April  1971).  (CONFIDENTIAL) 
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beam  active  optical  detection  systems  where  transmitter  and  receiver  are 
approximately  collocated.  However,  multiple  scattering  effects  are 
ignored  in  deriving  equation  (2). 

To  see  the  usefulness  of  equation  (2)  more  fully,  suppose  that  a 
pencil  beam  system  satisfying  the  general  requirements  of  the  previous 
paragraph  is  operating  in  an  aerosol  environment  and  that  it  is  desired 
to  calculate  the  aerosol  return  signal.  Assuming  that  the  aerosol 
distribution  is  known  in  terms  of  its  extinction  and  backscatter  coef¬ 
ficients,  equation  (1)  can  be  used  to  calculate  the  function  C(x)  once 
the  location  of  the  transceiver  and  the  orientation  of  its  pencil  beam 
influence  pattern  relative  to  the  aerosol  distribution  are  specified. 
Knowledge  of  the  range  sensitivity  curve  and  the  instantaneous  trans¬ 
mitted  power  then  permits  the  calculation  of  the  desired  return  signal 
by  using  equation  (2).  C(x)  is  system  independent,  being  determined 
solely  by  the  aerosol  distribution  and  the  encounter  geometry.  For  this 
reason,  C(x)  is  referred  to  here  as  the  aerosol  signature  (for  a  given 
encounter) . 

The  aerosol  signature  concept  easily  extends  to  wide  angle  influence 
pattern  systems,  since  these  can  be  viewed  as  superpositions  of  many 
pencil  beam  systems;  however,  it  is  unlikely  that  multiple  scattering 
effects  can  be  ignored  altogether  in  the  wide  angle  case.  To  the  extent 
that  such  effects  can  be  ignored,  the  aerosol  signature  concept  provides 
a  solution  to  the  problem  being  considered,  namely,  how  backscatter  data 
obtained  with  the  GaAs  laser  probe  can  be  made  to  apply  to  AOF  systems 
with  a  different  influence  pattern  and  different  transmitter  pulse 
characteristics.  One  needs  to  extract  the  aerosol  signatures,  C, 
contained  in  the  measured  aerosol  return  signals,  V,  these  being  related 
according  to  equation  (2),  with  K,  P,  and  R  understood  to  refer  to  the 
GaAs  laser  probe.  The  investigations  to  be  discussed  in  this  report  are 
concerned  mainly  with  how  to  accomplish  the  signature  extraction. 

There  is  a  simple  approximate  way  to  determine  the  aerosol  signature 
from  the  measured  return  signal.  The  laser  probe  features  short  trans¬ 
mitter  pulses  (5  to  11  ns  FWHM)  and  is  sensitive  to  aerosol  over  a  range 
interval  (in  excess  of  10  m)  consistent  with  fuzing  applications.  If 
one  approximates  the  transmitter  pulse  by  a  temporal  6-function,  then 
equation  (2)  shows  that  the  return  signal  is  proportional  to  the  product 
of  the  aerosol  signature  and  the  probe's  range  response  function.  Thus 
the  aerosol  signature  can  be  obtainSd  approximately  by  dividing  the 
return  signal  values  (expressed  as  a  function  of  range  using  the  trans¬ 
formation  t  2x/c)  by  the  corresponding  values  of  the  range  response 
function,  appropriately  scaled.  This  approximate  determination  sacri¬ 
fices  resolution  in  the  resulting  C(x).  Since  the  spatial  width  of  the 
transmitter  pulse  is  about  2  m,  the  C(x)  obtained  will  be  smeared  rela¬ 
tive  to  the  actual  C(x)  by  a  spatial  average  over  about  a  2-m  interval. 
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The  general  signature  extraction  problem  is  basically  to  solve 
equation  (2)  as  an  integral  equation  for  C(x),  given  V(t),  P(t),  and 
R(x),  and  to  analyze  the  errors  produced  in  the  solution  by  the  noise 
accompanying  V(t).  In  section  2,  the  problem  is  formulated  as  one  of 
deconvolution,  that  is,  inverting  a  convolution  operator,  and  a  useful 
point  of  view  for  analyzing  noise  effects  is  indicated.  In  addition, 
the  transmitter  pulses  and  range-response  function  of  the  GaAs  probe  are 
described,  and  a  number  of  idealized  aerosol  signatures  are  computed  and 
displayed  graphically. 

In  section  3,  a  formal  solution  of  the  deconvolution  problem  is 
given  in  terms  of  Fourier  transforms,  and  a  formulation  of  the  noise 
analysis  problem  is  developed  on  the  basis  of  the  formal  solution.  Then 
a  method  is  presented  for  the  numerical  solution  of  equation  (2)  invol¬ 
ving  the  use  of  the  discrete  Z-transform.  This  method  leads  to  a  highly 
efficient  and  accurate  solution  algorithm.  A  sample  computation  illus¬ 
trating  the  accuracy  is  given. 

Section  4  completely  analyzes  the  noise  induced  errors  in  the 
general  signature  extraction  process.  As  might  be  expected,  the 
critical  factors  for  such  errors  are  the  width  and  the  shape  of  the 
transmitter  pulse.  By  making  reasonable  assumptions  about  the  statis¬ 
tics  of  the  noise  present  in  the  backscatter  data,  the  errors  in  deter¬ 
mining  C(x)  by  deconvolution  can  be  expressed  in  the  form  of  a  signal- 
to-noise  ratio  (SNR).  The  SNR  of  C(x)  is  calculated  for  a  transmitter 
pulse  that  much  of  the  backscatter  data  were  obtained  with;  the  pulse  is 
somewhat  asymmetrical  (the  leading  edge  is  faster  than  the  trailing 
edge)  and  measures  7  ns  FWHM.  This  calculation  was  performed  by 
assuming  that  the  backscatter  data  were  bandpass  filtered  with  various 
upper  frequency  cutoffs  so  that  the  effect  of  data  smoothing  could  be 
studied.  The  results  show  that  the  data  in  question  cannot  be  profit¬ 
ably  treated  with  the  deconvolution  method  because  the  SNR's  of  C(x)  are 
unacceptably  low  for  reasonable  amounts  of  data  smoothing  and  the 
typical  SNR  of  the  backscatter  data  (about  15:1,  peak  signal  to  rms 
noise)  . 

The  foregoing  negative  conclusions  do  not  apply  to  backscatter  data 
that  are  now  being  obtained  (and  that  were  unavailable  when  this  study 
was  conducted) .  The  current  data  are  being  obtained  with  a  shorter  (5 
ns  FWHM) ,  more  nearly  symmetrical  transmitter  pulse  that  has  consider¬ 
ably  more  peak  power  than  was  previously  available.  The  combination  of 
increased  SNR  for  the  backscatter  data  and  changed  temporal  character¬ 
istics  of  the  transmitter  pulse  could  make  the  deconvolution  method 
useful  for  the  newer  data.  Some  speculation  along  these  lines  is 
offered  in  section  6. 


A  preliminary  discussion  is  given  in  section  5  of  an  alternative 
strategy  for  accomplishing  signature  extraction.  This  strategy  involves 
the  use  of  a  priori  knowledge  about  the  signature  and  parameter  estima¬ 
tion  techniques . 


2.  PRELIMINARY  DISCUSSION  OF  SIGNATURE  EXTRACTION  PROBLEM 

Before  considering  methods  for  signature  extraction,  it  is  useful  to 
express  the  ideas  surrounding  equation  (2)  in  the  language  of  linear 
system  theory  and  to  describe  in  more  detail  the  several  functions — 
P(t) ,  R(x)  ,  and  C(x) — that  enter  into  the  problem. 

2. 1  Linear  System  Description 

Equation  (2)  can  be  seen  to  express  V(t)  as  a  convolution 
integral  if  the  convention  is  adopted  of  putting  C  and  R  equal  to  zero 
for  negative  values  of  their  arguments.  Then  the  lower  limit  on  the 
integral  can  be  changed  to  -*  without  affecting  the  validity  of  the 
equation.  Define  a  new  function  h  by 

h(t)  =  K0( ct/2 ) R( ct/2 )  .  (3) 

Then  V(t)  can  be  expressed  as 

V(t)  =  P(t-T)h(T)  dx  ,  (4) 

or,  using  *  to  denote  the  convolution  operation, 

V  =  P  *  h  .  (5) 

Equations  (4)  and  (5)  can  be  interpreted  as  saying  that  V(t)  is  the 
output  of  a  linear  system  whose  input  is  P(t)  and  whose  impulse  response 
is  h(t).  Since  the  problem  of  determining  the  aerosol  signature  is 
essentially  that  of  determining  h,  it  may  be  said  that  our  problem  is  to 
deconvolve  V  with  respect  to  P.  Other  ways  of  phrasing  the  problem  are 
(  1  )  determine  the  impulse  response  of  the  system  from  knowledge  of  its 
input  and  output  and  (2)  de  . ermine  the  input  h(t)  that  produces  a  known 
output  V(t)  when  the  impulse  response  P(t)  is  known.  The  latter 
phrasing  follows  from  the  fact  that  P  *  h  =  h  *  P  and  is  a  convenient 
way  of  looking  at  the  problem  for  noise  analysis. 


2.2  Typical  Transmitter  Pulses,  Range  Law,  and  Signatures 


The  laser  probe  has  undergone  several  modifications  in  the 
course  of  its  history  to  improve  transmitter  pulse  and  receiver  noise 
characteristics.  Figure  1  shows  the  shape  of  a  typical  transmitted 
pulse  for  the  probe's  original  constitution,  measured  by  using  an  ap¬ 
proximately  100-ps  response  time  photodiode.  The  measurement  was  made 
by  intercepting  the  transmitted  pencil  beam  with  the  sensitive  photo¬ 
diode  surface  (an  Instrument  Technology  Ltd.  HSD-SO)  and  recording  the 
temporal  variation  of  the  induced  photocurrent.  The  peak  transmitted 
power  is  roughly  5  to  10  W,  depending  on  whether  a  polarizer  is  inserted 
in  the  beam;  the  latter  is  often  the  case  in  experimental  flight  tests 
because  backscatter  depolarization  effects  also  are  measured.  The  pulse 
is  about  11  ns  wide  at  the  half-maximum  points  and  rises  noticeably  more 
rapidly  than  it  falls.  The  pulse  shape  needed  in  equation  (4)  is, 
however,  the  result  of  filtering  that  which  is  shown  in  figure  1  with 
the  bandpass  characteristic  of  the  probe's  receiver,  which  has  a  band¬ 
width  of  approximately  200  MHz.  This  follows  because  the  signal,  V(t), 
is  thusly  band  limited  and  from  the  basic  algebraic  properties  of  the 
convolution  operation.  Little  change  of  pulse  shape  would  occur  in 
figure  1  due  to  such  filtering,  though. 

The  first  probe  modification  resulted  in  the  faster  pulse  shown 
in  figure  2.  The  pulse  shape  was  measured  by  reflecting  the  transmitter 
pulse  into  the  probe's  receiver  system,  using  a  properly  oriented  flat 
reflecting  surface,  and  recording  the  resulting  waveform  at  the  receiver 
output.  Thus  the  pulse  shape  shown  in  figure  2,  which  is  about  7  ns 
FWHM,  is  that  appropriate  for  use  in  equation  (4).  The  transmitted 
optical  pulse  actually  has  a  much  faster  leading  edge  than  shown  in 
figure  2.  Measurements  using  the  ' JD-50  photodiode  showed  a  rise  time 
of  about  100  ps,  which  is  the  photodiode  response  time.  These  measure¬ 
ments  and  the  technical  details  of  the  pulser  design  are  discussed  by 
Vanderwall  et  al.6  The  receiver  bandpass  of  the  probe  slows  this 
quickly  rising  pulse  to  the  extent  indicated  by  figure  2. 


6 Jonathan  Vanderwall ,  Walter  V.  Hattery,  and  Zoltan  G.  Sztankay , 
Subnanosecond  Rise  Time  Injection  Laser  Pulses,  Harry  Diamond 
Laboratories  HDL-TR-1697  (March  1975). 
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Figure  1.  Shape  of  typical  optical  pulse  from  GaAs  laser  pulser 
measured  with  100-ps  response  time  photodiode.  This  pulser  was 
original  transmitter  section  of  GaAs  laser  probe  for  aerosol 
backscatter  measurements . 
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Figure  2.  Shape  of  S^pSS  i- 

pulser.  Shape  was  ' ^tca^  receiver  essentially 

approximately  200-MHz  an  x  P  appropriately  oriented  flat 

collocated  With  transmitter  section 

of  Improv^SfcHaser  probe  for  aerosol  backscatter 


All  backscatter  data  available  at  the  time  of  this  investi¬ 
gation  were  obtained  with  the  pulse  shapes  of  figures  1  and  2.  Recent 
modifications  have  further  shortened  the  transmitter  pulses  (to  about  5 
ns  FWHM)  and  have  substantially  increased  the  power  output  level. 

The  range  response  of  the  probe  can  be  adjusted  by  varying  the 
angle  between  the  transmitter  pencil  and  the  receiver  field  of  view  and 
by  varying  the  separation  between  the  transmitter  and  the  receiver. 
Figure  3  shows  a  typical  range-response  curve,  one  that  has  been  used  in 
several  flight  tests.  Apart  from  its  precise  shape,  the  curve's  main 
adjustable  features  are  the  range  at  which  peak  response  is  obtained  and 
the  range  interval  over  which  the  system  is  sensitive. 


Figure  3.  Measured  range-response  curve,  R(x) ,  of  pencil  beam  GaAs 
laser  probe  for  aerosol  backscatter  measurements.  Curve  was  obtained 
by  measuring  peak  receiver  response  to  flat  reflecting  target,  oriented 
at  right  angles  to  and  filling  influence  pattern,  as  function  of  range 
from  transceiver  to  target.  Circled  points  indicated  measurements. 
Scale  of  R  has  been  normalized  so  that  response  at  7  m  is  1/(7  m)2. 


The  general  features  of  the  aerosol  signature,  C(x) ,  can  be 
determined  by  considering  a  somewhat  idealized  aerosol  density 
profile.  Suppose  that  a(x)  and  p(x)  are  as  shown  in  figure  4.  Then  C 
vanishes  for  x  xQ  .  For  xQ  <  x  <_  x1  , 

/*  a(s)  ds  =  ( a/l )  /*  (s  -  xQ^  ds  =  (1/2)(a/2.)  (x  -  xQ^ 2  , 


so  that 


For  x  >  x. , 


C(x)  =  (p/£) 


(x  ■  xo) 


-(a/l)  ^c-x^ 


o  ( s )  ds 


=  a(s)  ds  +  Jx  o  ds  =  ( 1/2 )  (a/l  )l2  +  a  (x  -  xyj 


so  that 


-a l  ~2o(x-x  ) 
C(x)  =  pe  alLe  1 


If  distances  are  measured  in  units  of  l,  then  the  three  results  for  C(x) 
become 


0  ,  for  x  £  x„  , 


)  X  (x~xn\  ^ 

C(x)  =  (  u  -  xfl j  e  V  u/  ,  for  xQ  <  x  £  x1  , 


,  for  x  >  x, 


Thus,  with  range  measured  in  units  of  l,  the  only  parameter  upon  which 
C(x)/p  depends  is  al .  Plots  of  C( x)  according  to  equations  (8)  to  (10) 
for  various  values  of  a l  are  given  in  figure  5.  The  particular  value  of 
P  determines  the  vertical  scale,  and  the  particular  value  of  a  for  a 
given  curve  in  the  figure  determines  the  horizontal  scale  by  determining 
l.  Put  in  a  more  direct  way,  given  the  p,  a,  and  £  of  interest,  one 
chooses  the  curve  corresponding  to  the  product  of  a  and  l ;  then  p  deter¬ 
mines  the  appropriate  vertical  scale  to  use  and  l  determines  the  corre¬ 
sponding  horizontal  scale. 


*0  *1  =  *0  +  >’ 


x 


Figure  4.  Illustration  of  idealized  aerosol  density 
profile  along  pencil  beam  influence  pattern  of  active 
optical  detection  system.  Profiles  of  extinction 
coefficient,  a,  and  backscatter  coefficient,  u,  are 
indicated  as  function  of  range,  x,  from  transceiver. 


Figure  5.  Aerosol  signatures,  C(x),  calculated  from  defining 
equation  (1)  for  aerosol  density  profile  shown  in  figure  4  for 
various  values  of  dimensionless  parameter  ai  (noted  near 
signature  peaks) •  Range  is  measured  in  units  of  l,  thickness  of 
aerosol  buildup  region. 
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Two  points  should  be  made  about  figure  5.  First,  the  corners 
in  the  signatures  at  xQ  and  x  are  due  to  the  corners  in  the  corre¬ 
sponding  aerosol  density  profile.  Since  the  latter  never  occur  in 
reality,  the  former  are  always  rounded  off  in  real  cases.  Second,  the 
mode  of  presentation  is  inappropriate  for  signatures  corresponding  to 
the  limiting  case  l  =  0  because  then  the  entire  horizontal  scale  shrinks 
to  a  single  point;  alternatively,  measuring  range  in  units  of  i  =  0  is 
nonsense.  For  this  limiting  case,  the  signature  is  given  by 


C(x)  = 


0  ,  for  x  <  xQ  =  , 


(x“xl) 


,  for  x  X  x,  , 


which  is  sketched  in  figure  6. 

One  further  point  should  be  made  concerning  the  aerosol  density 
profile  used  in  the  foregoing  calculations.  The  aerosol  is  assumed  to 
extend  infinitely  far  in  the  positive  range  direction,  and  this  exten¬ 
sion  never  occurs  in  reality.  If  the  range  at  which  the  aerosol  density 
begins  to  drop  out  is  denoted  by  x2> then  equations  (8)  to  (10)  remain 
valid,  except  that  in  equation  (10)  the  validity  is  restricted  to 

xx  <  X  <  x2  . 

The  net  effect  on  the  signatures  of  figure  5  is  that,  from  x2  to  the 
range  where  the  aerosol  finally  ceases  to  exist,  C(x)  falls  to  zero. 


j-Zoli-io) 


Figure  6.  Aerosol  signature  for  semi-infinite 
uniform  aerosol  distribution  with  abrupt 
buildup  at  distance  x  from  transceiver. 


3.  FORMAL  AND  NUMERICAL  SOLUTIONS  OF  DECONVOLUTION  PROBLEM 


3.1  Formal  Solution  and  Effect  of  Noise 


The  deconvolution  problem  can  be  solved  formally  in  a  straight¬ 
forward  manner,  starting  from  equation  (5).  The  Ftourier  transform  of 
this  equation  is 

V<f)  =  P(f)h(f)  ,  ( 13) 

where  f  is  the  frequency  variable  and  hats  are  used  to  denote  Fourier 
transforms.  Ttius 


h(f)  =  V(f)/P(f)  ,  (14) 

provided  that  P(f)  *  0.  The  time  functions  V(t),  P(t),  and  h(t)  each 
vanish  outside  a  certain  time  interval.  Their  vanishing  guarantees  the 
existence  and  the  continuity  of  their  Fourier  transforms.  If  f  is  a 
frequency  for  which  P(fQ)  =  0,  then  equation  (13)  shows  °  that 
V(fo)  =  0  as  well.  Consequently,  the  right-hand  side  of  equation  (14) 
takes  the  indeterminate  form  0/0  for  such  frequencies.  It  is 
nevertheless  true  that  h(fQ)  is  well  defined,  and  the  continuity  of  all 
the  transforms  involved  shows  that 

A  A  A 

h(f  )  =  lim  V(f)/P(f>  .  (15) 

o 


< 


It  is  well  known  that  time  functions  that  vanish  outside  some  time 
interval  give  rise  to  Fourier  transforms  that  can  vanish  at  isolated 
frequences  only.  There  is  therefore  no  possibility  in  equation  (15) 
that  P  will  vanish  in  a  neighborhood  of  fQ.  Thus,  if  equation  (14)  is 
understood  to  mean  the  limit  in  equation  (15)  at  those  isolated  frequen¬ 
cies  where  it  is  singular,  the  impulse  response  h(  t)  can  be  obtained 
from  it  by  using  the  inverse  Fourier  transform,  namely. 


h(t) 


(16) 


Equation  (16)  expresses  a  formal  solution  to  the  deconvolution 
problem.  To  implement  such  a  solution  numerically  involves  several 
difficulties,  not  the  least  of  which  is  handling  the  singular 
frequencies  at  which  V  and  P  vanish.  Any  numerical  algorithm  that 
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computes  the  values  of  V  and  P  from  data  giving  V(t)  and  P(t)  poten¬ 
tially  commits  large  errors  in  determining  the  ratio  V/P  near  singular 
frequencies.  These  errors  may  in  turn  contribute  significant  error  in 
the  numerical  evaluation  of  the  integral  in  equation  (16). 

A  related  difficulty  of  critical  importance  pertains  to  the 
effect  of  noise  on  the  signal  V(t)  in  producing  errors  in  the  deter¬ 
mination  of  h(t).  At  singular  frequencies  of  P,  the  Fourier  transform 
of  a  noisy  signal  V(t)  +  n(  t)  is  that  of  the  noise,  which  can  only 
fortuitously  turn  out  to  be  zero.  Near  such  frequencies,  the  integrand 
in  equation  (16)  actually  tends  to  infinity  if  V  is  understood  to  be 
that  of  the  noisy  signal .  Thus ,  the  potential  for  very  large  noise 
related  errors  in  determining  h(t)  is  clear. 

When  the  noise  that  accompanies  the  signal  V(t)  is  known  to  be 
a  stationary  and  Gaussian  process,  which  is  a  reasonably  good  assumption 
for  the  laser  probe,*  a  method  for  analyzing  the  resulting  error  present 
in  the  determination  of  h(t)  can  be  based  on  the  linear  system  interpre¬ 
tation  of  equation  (5). 

The  method  is  suggested  by  writing  equation  (5)  in  the  equiv¬ 
alent  form 

V  =  h  *  P  ,  (17) 

which  can  be  interpreted  as  saying  that  h  is  the  input  to  a  linear 
system  with  impulse  response  P  and  output  V.  Given  that  the  output  is 
accompanied  by  a  well-characterized  noise  process,  the  following  can  be 
asked:  What  noise  process  must  accompany  the  input,  h,  in  order  that 

the  given  output  noise  will  be  produced,  assuming  that  the  system,  P,  is 
internally  noiseless?  This  question  can  be  answered  in  explicit  mathe¬ 
matical  terms  when  the  noise  processes  involved  are  known  to  be  station¬ 
ary  and  Gaussian.  Let  <N>  and  S^f)  denote  the  mean  and  the  spectral 
density  of  the  desired  input  noise  process.  Then  the  mean  <n>  and  the 
spectral  density  sn(f)  of  the  output  noise  are  given  by 


<n>  =  <N>P (0 )  , 

s  (f)  =  |P(f) | 2 S  ( f )  .  (18) 

n  N 


* Virtually  all  the  significant  noise  sources  present  in  direct 
optical  detection  systems  are  of  the  Johnson  or  shot-noise  type,  except 
possibly  the  avalanche  multiplication  noise  that  arises  when  avalanche 
photodiodes  are  used  for  detection. 


It  follows  that 


S  (f)  =  s  (f)/|P(f)  I2  M9) 

N  n 

for  all  frequencies  such  that  P(f)  *  0.  The  singular  frequencies,  fQ, 
at  which  P  vanishes  are  isolated;  however,  they  may  exist,  and  unless 
sn(f)  tends  to  zero  sufficiently  rapidly  as  f  ♦  fQ,  SN(f)  must  become 
infinite  there. 

The  quantity  of  interest  is  actually 


/  Vf)  df  ' 

— <X> 


since  this  is  the  mean-squared  noise  level  associated  with  h;  denote 
this  quantity  by  <N2>.  If  <N2>  can  be  calculated,  it  would  provide  a 
precise  characterization  of  the  error  in  h(t)  produced  by  the  noise 
accompanying  the  signal  V(t).  From  this  standpoint,  infinite  singular¬ 
ities  in  ^(f)  are  acceptable  as  long  as  they  do  not  cause  the  integral 
giving  <N2>  to  diverge. 

To  progress  further  with  the  analysis,  some  way  of  dealing  with 
the  infinite  limits  on  the  integral  in  question  is  needed.  To  see  why, 
note  first  that  both  the  noise  spectrum  sn(f)  and  P(f)  refer  to  those 
that  are  seen  at  the  receiver  amplifier  output.  In  the  case 
of  P(f),  the  P(t)  from  which  P(f)  is  to  be  determined  is  the  result  of 
filtering  the  actual  transmitted  pulse  with  the  receiver  bandpass  char¬ 
acteristic.  .  Since  it  is  clear  that  P(f)  +  0  as  |f|  ■*  °°,  it  follows 
that  |P(f)|~  -►  00  in  the  same  limit.  Thus,  for  the  integral  in  question 
to  converge,  it  is  necessary  that  sn(f)  tend  to  zero  sufficiently 
rapidly  as  |f|  *  °°.  Since  no  measurements  of  sn(f)  or  |P(f) |  are  going 
to  reveal  the  precise  analytic  behavior  of  these  quantities  for 
large  |f|,  recourse  to  some  other  procedure  that  fixes  the  asymptotic 
behavior  of  the  integrand  is  necessary.  A  reasonable  procedure, 
although  certainly  not  the  only  one,  is  to  cut  off  the  integration  at 
some  frequency,  fc,  that  is,  to  consider  only 
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This  is  tantamount  to  inserting  an  ideal  low-pass  filter  between  the 
amplifier  output  and  our  observation  point,  which  procedure  is  justifi¬ 
able  if  it  can  be  argued  that  the  chosen  cutoff  frequency  introduces 
negligible  distortion  in  the  return  pulse  from  the  aerosol. 

The  foregoing  analysis  is  pursued  in  section  4,  where  analyt¬ 
ical  estimates  of  the  effect  of  noise  are  obtained.  We  now  turn  to  the 
numerical  solution  of  the  deconvolution  problem. 

3.2  Numerical  Deconvolution 

The  inversion  of  the  convolution  operator  by  Fourier  trans¬ 
forms,  while  convenient  for  analytical  purposes,  does  not  lead  to  the 
most  efficient  numerical  solution  of  the  problem.  To  numerically  imple¬ 
ment  such  a  solution  most  efficiently  would  involve  the  use  of  the  fast 
Fourier  transform.  Another  approach,  which  is  ideally  suited  to  the 
digital  form  of  the  return  pulse  data,  uses  Z-transform  techniques.  The 
Z-transform  is  a  discrete  counterpart  of  the  Laplace  transform  and  has 
numerous  applications  to  digital  systems  and  sampled  data.7  To  apply  Z- 
transform  methods  to  the  deconvolution  problem  requires  first  that  the 
problem  be  cast  in  discrete  form.  The  method  leads  to  a  highly  effi¬ 
cient  numerical  solution. 

3.2.1  Discrete  Formulation 

The  transmitted  pulse,  P(t)  ,  and  the  return  signal,  V(t),  are 
observed  for  N  periodically  sampled  times,  A  seconds  apart,  in  the 
interval  [o,t],  where  (N  -  1)A  =  T.*  In  this  case,  equation  (41  is 
equivalent  to 

,t 

V(t)  =  /  PfT)h(t-T)  dx  ,  (20) 

0 


where  the  commutativity  of  convolution  has  been  used  and  it  is  assumed 
that  P(t)  =  0  for  t  <  0.  For  the  rest  of  the  formulation,  it  is  neces¬ 
sary  to  assume  also  that  P(0)  t  0;  this  restriction  may  always  be 
achieved  by  finite  time  translation.  A  convenient  quantization  for  P 
and  h  is 

P(t)  =  P(nA)  ,  for  nA  <  t  £  (n  +  1)A  ,  (21) 

/E.  I .  Jury,  Theory  and  Application  of  the  Z-transform  Method ,  John 
Wiley  &  Sons,  Inc,,  New  York  (1964). 

*The  digitization  of  the  data  results  in  N  =  512 ,  with  A  typically 
about  200  ps. 
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and 
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Since  P(0)  was  assumed  nonzero,  it  follows  that  the  determinant  of  P 
does  not  vanish  so  that  the  matrix  P  has  a  unique  inverse  P-1 .  The 
solution  of  equation  (27)  is  therefore 


h  = 


(28) 


3.2,2  Z-Transform  Deconvolution 


The  foregoing  discrete  formulation  shows  that  the  numerical 
deconvolution  problem  is  essentially  to  find  the  inverse  of  a  large 
matrix  or,  equivalently,  to  solve  a  large  number  of  simultaneous  linear 
equations,  namely,  equation  (27).  The  Z-transform  method  can  be  applied 
to  equation  (27)  and  results  in  a  rather  high-speed  inversion  of  P.  To 
illustrate  the  results  obtainable  in  this  manner,  we  show  in  figures  7 
to  10  the  results  of  a  sample  calculation. 

.  Equations  (6)  and  (7)  were  used  with  the  parameter  values  a  = 
0.2  m  ,  U  =  0.01  m  sr  ,  Xy  =  2.75  m,  and  x^  =  5.25  m  to  compute  the 
aerosol  signature  shown  in  figure  7.  The  parameter  values  correspond  to 
a  potential  encounter  with  a  dense  water  cloud.  The  equation 


{0  ,  x  <  2.67  m  , 

1.83(1  -  2.67/x) ( 1/x2 )  ,  2.67  m  <  x  <  5.89  m  ,  (29) 

1/x2  ,  x  >  5.89  m  , 


was  used  to  compute  the  range- response  function  R(x)  ,  which  fits  the 
measured  range-response  function  of  figure  3  with  fair  to  good 
accuracy.  The  product  of  R(x)  and  the  aerosol  signature  (which  is 
proportional  to  h)  is  shown  in  figure  8.  A  unit  amplitude  transmitter 
pulse,  u(t),  having  the  form 


I 


( 1/4) [ 1  +  cos  (vt  -  n)]2 


0  <  t  <  2n/v  , 


u(  t) 


0  ,  otherwise 


was  then  assumed;  the  frequency  parameter  v  (»  471  MHz)  was  set  by  the 
condition  ctt/v  =  2  m  (cir/v  is  the  spatial  basewidth  of  u(  t) ) .  This 
model  pulse  roughly  approximates  the  measured  transmitter  pulse  of 
figure  2,  apart  from  the  origin  of  the  time  axis,  and  was  used  to  calcu¬ 
late  the  return  signal  shown  in  figure  9  from  equation  (23).  In  equa¬ 
tion  (23),  h  was  taken  as  given  by  the  CR  product  of  figure  8;  that  is, 
the  proportionality  constant  K  was  set  equal  to  1,  and  A  was  set  at  2/3 
ns.  Finally,  the  Z-transform  method  was  used  to  deconvolve  the  return 
signal  of  figure  9,  that  is,  to  solve  equation  (27)  for  h  =  CR,  and  the 
result  is  shown  in  figure  10. 

The  agreement  of  figures  8  and  10  is  excellent,  as  can  be 
verified  by  superimposing  tracings  of  the  two  waveforms.  The  relative 
shift  in  range  between  the  waveforms  is  not  a  real  discrepancy;  the 
shift  is  caused  by  the  need  for  shifting  the  transmitter  pulse  in  time 
so  that  P(0)  t  0,  a  condition  required  for  the  deconvolution  compu¬ 
tations.  The  agreement  obtained  shows  that  noiseless  signals  can  be 
deconvolved  numerically  with  considerable  accuracy.  Moreover,  the 
machine  calculations  needed  to  produce  figure  10  were  sufficiently  rapid 
to  make  the  deconvolution  of  large  data  banks  of  aerosol  return  signals 
feasible.  Unfortunately,  the  effect  of  adding  even  very  small  amounts 
of  simulated  noise  to  the  return  signal  is  catastrophic.  To  the  issues 
surrounding  this  fact  we  now  turn. 
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Figure  8.  Product  of  aerosol  signature  of  figure  7 
and  range-response  function,  R(x),  computed  from 
equation  ( 29) . 


Figure  9.  Computed  aerosol  return  signal,  V, 
according  to  equation  (23)  for  h  =  CR(K=1)  as  given 
in  figure  8,  A  =  2/3  ns,  and  P(t)  =  u(t)  =  (t/4)[1  + 
cos  (vt-n)]2  (0  <_  t  <_  2n/v);  v  (»  471  MHz)  is  set 
by  cfl/v  =  2  m. 
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Figure  10.  Product  of  aerosol  signature  and  range 
response  function  as  determined  by  Z-transform 
deconvolution  of  return  signal  of  figure  9, 


4.  EFFECT  OF  NOISE  IN  SIGNATURE  EXTRACTION 


4.1  Numerical  Effect 


In  the  measurements  of  aerosol  return  signals  using  the  GaAs 
laser  probe,  various  noise  sources  combine  to  produce  a  typical  SNR 
(peak  signal  divided  by  rms  noise)  for  all  but  the  most  recently 
obtained  data  of  about  15.  To  determine  the  effect  of  such  a  noise 
level  on  the  sample  numerical  deconvolution  just  discussed,  the  return 
signal  of  figure  9  was  corrupted  with  simulated  noise,  giving  an  overall 
SNR  of  15,  and  subjected  to  the  Z-transform  deconvolution  algorithm. 
The  resulting  deconvolved  signal  is  shown  in  figure  11.  The  effect  of 
the  noise  was  catastrophic.  The  simulated  noise  was  produced  by 
filtering  random  numbers  with  a  third-order  Butterworth  filter  with  an 
upper  frequency  cutoff  of  200  MHz  (the  approximate  bandwidth  of  the 
amplifiers  used  in  the  laser  probe’s  receivers)  and  then  scaling  so  that 
the  rms  value  was  1/15  the  peak  value  of  the  return  signal.  To  deter¬ 
mine  if  the  catastrophe  would  disappear  for  higher  SNR's,  the  decon¬ 
volution  computation  was  repeated  for  an  SNR  of  150.  The  result  was 
similar  to  that  in  figure  11. 
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Figure  11.  Result  of  deconvolving  return  signal  of 
figure  9  after  it  has  been  corrupted  with  simulated 
noise  giving  overall  signal-to-noise  ratio  of  15. 
Noise  bandwidth  is  approximately  200  MHz. 


4. 2  Analytical  Estimation  of  Noise  Effects 


This  section  gives  an  analysis  of  the  effect  of  noise  by 
using  the  ideas  and  the  assumptions  of  section  3.1.  The  result  of  the 
analysis  is  that  a  good  understanding  of  the  phenomenon  shown  by  figure 
11  is  obtained,  and  the  means  for  overcoming  the  noise  problem  are 
indicated. 


The  starting  point  is  the  equation 


<N2> 


8n(f)/|P(f) 


df 


(30) 


for  the  mean-squared  noise  level  associated  with  h,*  which  is  propor¬ 
tional  to  the  product  of  the  aerosol  signature  and  the  range- response 
function.  Equation  (30)  follows  from  the  basic  relation 


<N2> 


II  s»|f|  a£  - 


(31) 


equation  ( 19) ,  and  the  method  chosen  to  fix  the  asymptotic  behavior  of 
the  integrand,  namely,  to  cut  off  the  integral  at  an  unspecified 
frequency  fc.  Although  fc  is  regarded  as  a  variable  parameter  in  what 
follows,  it  cannot  be  chosen  too  small  since  it  must  correspond  to 
enough  bandwidth  to  pass  the  return  signal  with  good  fidelity. 


To  make  an  explicit  evaluation  of  equation  (30)  possible,  it 
is  assumed  that  sn(f)  is  constant  over  the  frequency  interval  -f  <  f  j<_ 
fc  an<^  that  the  P(t)  from  which  P(f)  is  determined  is  given  by  the 
analytical  model 

^PQ  cos2  ( 7i  t/  2T )  ,  -T  <_  t  <_  T  , 


P(t)  =< 


(32) 


^0  ,  |t |  >  T  . 


* Equation  (30)  actually  gives  the  mean-squared  noise  level 
associated  with  the  result  of  convolving  h  with  the  impulse  response  of 
a  filter  with  sharp  lower-  and  upper- frequency  cutoffs  at  -fQ  and  fc, 
respectively . 
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A  plot  of  the  model  pulse  superimposed  on  the  measured  transmitter  pulse 
of  figure  2  is  given  in  figure  12.  PQ  and  T  were  chosen  to  give  the 
model  pulse  the  same  peak  value  and  the  same  FWHM  as  the  measured 
pulse.  The  model  pulse  is  only  a  rough  approximation  to  the  measured 
one  primarily  because  the  measured  pulse  is  not  symmetric  about  its 
peak.  The  model  pulse  has  the  advantage  of  being  easily  transformed 
mathematically  in  the  manner  required  for  the  calculation  of  <N2>, 
however . 
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Figure  12.  Comparison  of  shape  of  model  transmitter 
pulse  P(t)  =  PQ  cos2  (irt/2T)  with  that  of  measured 
transmitter  pulse  of  figure  2;  T  was  chosen  to  give 
model  pulse  same  FWHM  as  measured  pulse,  and  both 
pulses  are  aligned  to  agree  at  their  peaks. 
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The  assumption  that  sn( f )  is  constant  over  the  frequency 
interval  -f  <_  f  _<  fc  implies  that  the  mean-squared  noise  level  <n2>  on 
the  return  signal  is  related  to  sn  by 

sn  =  <nz>/ 2fc  ;  (33) 

the  return  signal  being  referred  to  is  the  actual  return  signal  after 
being  put  through  an  ideal  low-pass  filter  with  cutoff  frequency  fc. 
The  Fourier  transform  of  the  model  transmitter  pulse  is 


P(f)  = 


sin  2rrfT 


1  -  ( 2f T) ‘ 


so  that 


<n2>  fctrfT  El  -  (2fT)2])2 


|P(f)  I2  2f  P2T2 
c  o 


sin  2itfT 


By  changing  the  independent  variable  to  x  =  2fT,  <N2>  can  be 


written  as 


2f  P2T3 
c  o 


rx)2  (l  -  x2)  2 


sin‘  itx 


Denote  the  integrand  by  I(x)  and  2fcT  by  xc.  Then 


<N2>  = 


<n2> 

2f  P^T3  Jo 

c  o 


I(x)  dx 


I(x)  is  everywhere  positive  and  has  infinite  singularities  at  x  =  2,  3, 
.  .  .  ;  it  is  readily  verified  that,  although  I  takes  the  indeterminate 

form  0/0  when  x  =  0  or  1 ,  I(0)  =  1  and  1(1)  =  4.  A  plot  of  I(x)  is 
given  in  figure  13;  the  x-axis  has  also  been  labeled  with  the  corre¬ 
sponding  frequency  variable,  assuming  that  T  =  7  ns.  It  is  straight¬ 
forward  to  verify  that 


lim  /  C  I(x 
x  +2  Jo 


)  dx  = 


Figure  13.  Sketch  of  graph  of  I(x)  =  ("xlMI  -  x^JVsin^  "X 
Frequency  axis  is  determined  by  x  =  2fT,  with  T  =  7  ns. 


Consequently,  for  the  model  being  used,  it  makes  no  sense  to  consider 
cutoff  frequencies  beyond  T-L  13  143  MHz.  For  xc  <  2,  define 


I(x)  dx  , 


(39) 


which  is  plotted  versus  xc  in  figure  14.  The  integral  was  evaluated 
numerically  by  using  the  trapezoidal  rule  and  an  x-axis  grid  spacing  of 
0.1.  (This  approximation  slightly  overestimates  the  integral.)  Since 


<N2> 


(V^)  J(xc) 


(40) 


(eq.  33  has  been  used  here),  the  curve  of  figure  14  shows  the  variation 
of  <N2>  with  the  cutoff  frequency. 


The  mathematical  origin  of  the  divergent  behavior  of  <N2> 
as  fc  +  T-  is  the  vanishing  of  P( f )  as  f  +  T*  .  Any  transmitter  pulse 
whose  Fourier  transform  has  real  zeros  will  cause  the  type  of  divergence 
indicated  in  figure  14.  This  fact  can  be  seen  as  follows.  Let  f  denote 
the  smallest  positive  zero  of  P(f).  Since  P(f)  is  an  even  function,  it 
follows  that 


<N‘ 


>->f: 


s  (f)/|P(f) I2  df  =  g 
n 


('•) 


(41  ) 


Since  P(f)  is  necessarily  an  analytic  function  for  all  f,  it  has  the 
Taylor  expansion 


_(f-f)2+...,  (42) 

f=f 

which  converges  for  all  f.  In  a  sufficiently  small  neighborhood 
of  f,  P(f)  can  be  uniformly  approximated  by  an  expression  of  the 
form  a(f  -  f)  ,  where  "a"  is  a  nonzero  constant  and  m  >_  1  is  an 
integer.  Thus  if  sn(f)  is  merely  finite  and  nonzero  in  a  neighborhood 
of  f,  it  follows  that 

lim  g  ^f  =  €0  •  (  43 ) 

f  f 
c 

Equations  (41)  and  (43)  show  that  <N2>  00  as  f  f . 


P(f)  =|| 


-  1  d2P 

-  (f  -  f)  +  -  - 

f=f  2  df2 
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The  measured  transmitter  pulse  shown  in  figure  2  has  been 
numerically  Fourier  transformed  for  frequencies  up  to  250  MHz,  and  the 
transform  has  been  found  to  have  no  zeros  in  this  range.  Consequently, 
the  divergence  problem  will  not  occur  for  this  transmitter  pulse. 


It  can  be  argued  that  the  model  being  used  is  approximately 
correct  for  all  values  of  f  <  T-^  .  The  validity  of  the  model  with  a 
restricted  range  for  fc  can  be  seen  by  considering  the  problem  in  a 
broader  perspective.  Let  h^  denote  the  impulse  response  of  an  arbitrary 
low-pass  filter.  The  mean-squared  noise  level  accompanying  h  *  h£  i r, 
given  by  an  expression  similar  to  that  of  equation  (30),  namely. 


s  (f) 
n 

|P(f)  I2 


Ih^f)  I2 


df  . 


This  result  follows  by  applying  the  arguments  surrounding  equations  (17) 
to  (19)  to  the  equation  V  *  h£  =  (h  *  h£ )  *  P  and  then  noting  that  the 
spectral  density  of  the  noise  accompanying  V  *  h£  is  given  in  terms  of 
that  for  V  by  |h^ ( f )  |2sn{ f ) .  By  choosing  an  ideal  low-pass  filter  with 
a  sharp  cutoff  at  fQ,  we  effectively  took 


h  ( t)  =  2f  (sin  2itf  t)/2rr  f  t  =  h  (t)  (44) 

x.  c  c  c  c 

since  |h  ( f ) | 2  =  1  for  |f|  £  f  and  vanishes  otherwise.  For  the  general 
filter,  the  basic  requirement  for  the  analysis  is  that  the  integral 
expression  preceding  equation  (44)  be  finite.  Divergence  can  come  from 
two  sources:  the  infinite  limits  on  the  integral  and  any  zeros 

of  P(f).  The  infinite  limits  are  handled  by  choosing  filters  with 
appropriate  asymptotic  behavior  as  f  -*•  ±®,  such  as  those  with  sharp 
cutoffs.  The  zeros  of  P(f)  can  be  handled  only  by  choosing  filters  such 
that  hg_  vanishes  sufficiently  fast  at  the  zeros  to  overcome  the  diver¬ 
gence  that  is  automatically  present  in  equation  (41)  as  fc  tends  to  a 
zero  of  P(f).  Therefore,  the  basic  requirement  for  the  ideal  low-pass 
filter  with  sharp  cutoff  frequency  fc  is  that  fQ  be  less  than  the  small¬ 
est  positive  zero  of  P(f).  The  last  point  is  important  for  under¬ 
standing  the  behavior  shown  in  figure  11. 

We  now  proceed  with  the  analysis  for  the  model  transmitter 
pulse,  keeping  the  points  just  raised  in  mind. 

Let  h^  denote  the  maximum  value  of  h  *  he •  Then  according  to 
equation  (40),  the  SNR  necessarily  connected  with  the  determination  of 
h  *  hc  is 


39 


(SNR)h*hc  =  ^P0T3/yfnj(xc|1/2 


(45) 


The  SNR  associated  with  V  *  he  is 


(SNRJv*^  «  Vil/(<n2>)1/2 

=  VJ/  (2fcsnj/2 


(46) 


where  V'  is  the  peak  value  of  V  *  he*  Since  the  transmitter  pulse  has 
very  short  duration,  a  large  error  will  not  result  from  using  the 
approximate  relationship  »  Eh^,  where  E  =  PQT  is  the  total  energy  in 
the  transmitted  optical  pulse  described  by  equation  (32).  Equations 
(45)  and  (46)  then  combine  to  yield 

( SNR )  h*  “  [xc/j(xcj]  1/2  ( SNR)  .  (47) 


If  the  SNR's  of  V  and  V  *  hc  are  compared,  one  finds  approximately  that 

( SNR) v#hc  -  (fa/fcy/2(SNR)v  ,  (48) 

where  f  is  the  receiver  amplifier  bandwidth.  Combining  equations  (47) 
and  (48)  results  in 

( SNR ) h* he  •  [xc/j£c)]1/2  (fa/fC)  1/2(SNR)v  ,  (49) 

which  is  the  main  result  of  the  analysis. 


Equation  (49)  indicates  how,  depending  on  the  chosen  cutoff 
frequency,  noise  errors  in  the  return  signal  V  are  translated  into  noise 
errors  in  h  *  hc*  The  factor  [xc/ J(xc) ]  versus  xc  is  plotted  in 
figure  15.  It  shows  the  progressive  degradation  in  (SNR)h*hc  compared 
with  (SNR)v*hc  as  the  cutoff  frequency  increases  toward  “143  MHz, 
where  (SNR)h*hc  vanishes.  The  situation  improves  significantly  if  the 
measured  transmitter  pulse  of  figure  2  is  used  for  the  analysis.  Figure 
16  plots  the  SNR  penalty  factor  of  this  case,  namely, 
{ SNR)  ft*  he/ ( SNR)  v*hc  •  Also  included  in  the  figure  is  a  segment  of  the 


40 


curve  of  figure  15  for  comparison.  The  SNR  penalty  factor  curve  for  the 
measured  transmitter  pulse  was  obtained  from  the  numerically  determined 
Fourier  transform  of  the  pulse  in  figure  2  by  evaluating  the  formula 


SNR  penalty  factor 


f 

c 


|P(0)  |2 


df 

|P(f) I 


2 


i/2 


(50) 


numerically  for  various  values  of  f  . 

If  f  a  f  a  200  MHz,  figure  16  indicates  that 
c  a 


(SNR)h*hc  »  0.15(SNR)V  . 


Thus  for  the  15-to-1  SNR  typical  of  much  of  the  aerosol  return  signal 
data,  the  SNR  of  h  *  hc  is  predicted  to  be  about  2.  If  the  aerosol 
return  signal  had  an  SNR  10  times  this  typical  value,  the  SNR 
of  h  *  hc  would  be  about  20.  These  predictions  apparently  contradict 
the  results  of  figure  11  and  the  corresponding  results  for  a  simulated 
return  signal  with  a  150-to-1  SNR.  There  are  important  differences 
between  the  transmitter  pulses  for  the  two  cases,  however. 

11  were  obtained  tor  a  transmitter 


1  11 /v  * 

(51) 


where  PQ  is  the  peak  transmitted  power  and  v  *  471  MHz.  By  a  somewhat 
tedious  calculation,  the  Fourier  transform  of  can  be  found  to  be 


pulse,  P^,  given  by 


The  results  of  figure 
given  by 

(po/4)(  1  +  cos  vt)2  ,  |t  | 


Pu(t)  = 


0  ,  1 1 1  >  u/v  , 


Pn<f> 


3P  IL  ^x 

o  v  irx 


f  4  (' 


(52) 


where  x  =  2»f/\>.  The  expression  for  P[  i  has  an  infinite  number  of 

zeros,  namely,  for  f  =  ±3^/2 k ,  ±4v/2",  ....  The  smallest  positive 

zero  is  f  =  3v/2u  ■  225  MHz. 


Figure  15.  Signal-to-noise-ratio  penalty  factor  (xc/J(xc)] 
versus  xc  and  f  for  model  transmitter  pulse  of  equation  (32). 
J(xc)  is  defined  by  equation  (39)  and  figure  13,  xc  =  2fcT,  and 
T  =  7  ns  for  fc  axis  labeling. 
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Figure  16.  Signal-to-noise  ratio  penalty  factor  versus  f  for  measured 
transmitter  pulse  of  figure  2  and  segment  of  signal-to-noise  ratio 
penalty  factor  curve  of  figure  15  for  comparison. 

The  Butterworth  filter  used  to  establish  the  frequency  band 
for  the  simulated  noise  that  produced  the  results  of  figure  11  does  not 
have  a  vanishing  frequency  response  at  225  MHz.  Therefore,  in  ac¬ 
cordance  with  previous  discussion,  divergent  behavior  in  the  noise 
accompanying  the  deconvolved  signal  is  expected,  regardless  of  the  SNR 
of  the  aerosol  return  signal. 

The  unacceptably  low  SNR  of  h  *  he  (■  2)  predicted  with 
equation  (50)  for  fc  =  200  MHz  and  a  measurement  SNR  of  15  to  1  can  be 
improved  by  lowering  fc.  This  tack  does  not  lead  to  acceptable  results, 
however,  because  reasonably  good  SNR's  are  obtained  only  for  values  of 
fc  that  significantly  affect  the  shape  of  the  return  pulses.  The  simple 
signature  extraction  procedure  of  assuming  that  the  transmitted  pulse  is 


approximately  a  temporal  fi-function  gives  an  approximation  of  h  that  is 
smeared  by  a  roughly  2-m  averaging  of  the  desired  impulse  response. 
This  approximation  is  corrupted  also  by  the  measurement  noise,  but  the 
SNR  is  the  same  for  h  as  for  the  measured  signal  V.  Thus  the  simple 
procedure  and  the  full  deconvolution  procedure  with  its  parameter  f  can 
be  compared  as  follows.  The  former  gives  an  SNR  equal  to  that  of  the 
return  pulse  data  and  a  signature  shape  distorted  by  a  roughly  2-m 
averaging.  In  the  latter,  both  the  SNR  and  the  signature  shape 
distortion  depend  on  fc  in  the  manner  of  a  tradeoff,  and  there  is  no 
advantage  in  using  the  more  complex  procedure  if  an  improvement  of  the 
2-m  averaging  and  an  acceptable  SNR  cannot  be  obtained  with  it.  Our 
results  show  that  such  improvement  and  an  acceptable  SNR  cannot  be 
obtained  for  the  15-to-1  measurement  SNR  typical  of  much  of  the  back- 
scatter  data.  Nonetheless,  a  sufficiently  high  measurement  SNR,  which 
may  now  be  available  due  to  recent  probe  modifications,  would  change 
this  conclusion. 

Several  additional  approaches  could  be  pursued  to  improve 
signature  extraction  accuracy.  One  could  attempt  to  modify  the  trans¬ 
mitter  pulse  shape  so  as  to  improve  the  SNR  penalty  characteristic.  One 
could  also  replace  the  sharply  cutoff  observation  filter  with  a  more 
general  type  and  seek  to  optimize  the  results  as  a  function  of  the 
filter  characteristic.  Certainly  the  most  direct  approach  is  to  further 
increase  the  peak  transmitted  power  so  as  to  improve  the  measurement 
SNR.  Finally,  the  use  of  a  priori  knowledge  about  C(x)  together  with 
parameter  estimation  techniques  could  be  pursued.  The  latter  approach, 
which  constitutes  a  new  strategy,  is  given  a  preliminary  discussion  in 
the  next  section. 

5.  ALTERNATIVE  APPROACH  THROUGH  PARAMETER  ESTIMATION 

The  problem  whose  discussion  has  occupied  the  bulk  of  this  report  is 
an  example  of  what  is  known  mathematically  as  an  ill-posed  problem,  that 
is,  one  whose  solution  does  not  depend  continuously  on  the  given  data. 
For  such  problems,  there  is  no  guarantee  that  reduction  of  the  errors  in 
the  given  data  will  reduce  the  error  in  the  solution,  and  this  lack  of 
guarantee  would  seem  to  indicate  that  the  practical  obtainment  of  solu¬ 
tions  to  ill-posed  problems  is  a  matter  of  fortuity.  Such  a  conclusion 
is  false,  however,  because  solution  methodology  is  available.  Ill-posed 
problems  are  currently  of  considerable  interest  in  a  number  of  applica¬ 
tion  areas,  particularly  in  the  geophysical  interpretation  of  seismic 
data.  The  general  approaches  to  a  practical  solution  include  the  selec¬ 
tive  reduction  of  the  information  sought  and  the  use  of  any  a  priori 
knowledge  concerning  the  object  under  study  to  provide  further  con- 
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straints  on  the  problem.  Our  use  of  a  low-pass  observation  filter  is  an 
example  of  the  selective  information-reduction  strategy.  This  section 
of  the  report  discusses  a  preliminary  investigation  of  applying  a  priori 
knowledge  together  with  the  use  of  parameter  estimation  techniques  as  an 
alternative  approach. 

The  main  source  of  a  priori  knowledge  for  the  signature  extraction 
problem  is  contained  in  equation  (1),  which  defines  the  signature  in 
terms  of  the  physical  properties  0  and  p  of  the  aerosol .  Let  us  model 
the  aerosol  distribution  in  the  influence  pattern  of  the  GaAs  probe  as  a 
sequence  of  layers  of  nonzero  thickness,  within  each  of  which  the 
physical  parameters  o  and  U  are  constant,  but  allowing  that  the  param¬ 
eter  values  in  each  of  the  layers  be  arbitrary.  This  distribution  of 
parameter  values  will  be  the  problem  solution  that  we  propose  to  deter¬ 
mine  from  measured  aerosol  return  signals.  How  can  this  determination 
be  made? 

Suppose  that  the  aerosol  distribution  model  has  M  layers,  so  that 
the  2M  parameters  c^,  u^,  0^,  .  .  .,  oM,  uM,  which  we  denote  vector- 

ially  by  o,  u,  are  sought.  By  using  equation  (1)  and  the  model,  the 
corresponding  aerosol  signature  C{x,g,M)  can  be  calculated  analytically 
and  depends  on  our  2M  parameters.  By  next  using  equation  (2),  a 
similarly  parameterized  return  signal  V(t,g,p)  can  be  calculated.  The 
final  step  is  to  seek  the  parameter  values  that  minimize,  in  some  sense, 
the  difference  between  V(t,g,g)  and  the  measured  return  signal. 

There  is  no  question  concerning  the  existence  of  solutions  to  the 
type  of  minimization  problem  outlined  since  the  dependence 
of  V(t,o,y)  on  o  and  y  is  continuous,  and  attention  can  be  restricted  to 
a  compact  region  of  the  2M-dimensional  parameter  space  in  seeking  solu¬ 
tions  (because  a  priori  we  know  reasonable  bounds  on  the  physical  pa¬ 
rameters)  .  Difficulties  could  arise,  however,  if  several  relative 
minima  are  present,  or  in  implementing  a  numerical  solution  for  fairly 
large  values  of  M. 

The  approach  just  outlined  makes  substantial  use  of  a  priori  infor¬ 
mation  about  C(x)  and,  if  successful,  will  yield  a  mathematical  formula 
for  a  continuous  solution  to  the  signature  extraction  problem.  More¬ 
over,  the  method  will  actually  give  approximate  extinction  and  back- 
scatter  coefficient  profiles  (along  the  influence  pattern)  so  that  more 
information  than  is  contained  in  the  aerosol  signature  may  be 
obtained.  The  validity  of  the  solution  will,  of  course,  be  limited  by 
the  noise  accompanying  the  measured  signal,  but  this  limitation  is  no 
worse  than  the  measurement  noise.  The  potential  lack  of  uniqueness  of 
the  solutions  could  pose  a  more  serious  problem  if  all  but  one  of  them 
could  not  be  ruled  out  on  a  priori  grounds'. 
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5.1 


Estimation  of  o  and  y  for  Uniform  Aerosol 


To  test  the  utility  of  the  foregoing  approach,  the  following 
problem  was  considered.  Suppose  that  the  GaAs  probe  is  fully  immersed 
in  a  uniform  density  aerosol  with  constant  values  of  a  and  y  through¬ 
out  .  Then 


_ ,  .  -2ox  ,  _  _ , 

C(x)  =  ye  .  (53) 

If  a  measured  return  signal  is  known  to  have  been  obtained  under  such 
circumstances,  then  the  problem  of  determining  the  values  of  a  and  y 
from  the  measured  signal  can  be  posed.  This  is  the  simplest  conceivable 
problem  of  the  type  being  considered.  We  discuss  its  solution  in  the 
discrete  framework  of  section  3.2.1. 

We  assume  that  the  transmitted  pulse  and  the  measured  return 
signal  are  observed_for  N  periodically  sampled  times,  A  seconds  apart, 
in  the  interval  [0,T],  where  T  =  (N  -  1)A.  The  quantization  chosen  by 
equation  (21)  gives  rise  to  the  N  x  n  matrix  P  of  equation  (26),  while 
equation  (24)  gives  the  corresponding  discrete  representation  of  the 
measured  return  signal.  Equation  (24)  defines  an  N-component  row  vector 
V'  so  that  the  corresponding  column  vector  is  denoted  by  V.  Let  h(a,y) 
denote  the  column  vector  with  components  ye  acn  R(cnA/2 ) ,  n  =  0,  1, 

•  •  .,  N-1.  Then  the  column  vector  representation  of  our  param¬ 

eterized  return  signal  V(t,a,y)  is  APh(a,y).  Thus  we  seek  to  minimize, 
in  some  sense,  the  vector  APh(a,y)  -  V. 

Tbe  minimization  problem  is  considered  for  a  general 
quadratic  cost  function 

J(o,y)  =  (APh(o,y)  -  V) ’A(APh(a,y)  -  V)  ,  (54) 

where  A  is  any  fixed  positive  definite  N  *  N  matrix;  the  prime  notation 
indicates  the  transpose  operation.  The  problem  is  to  find  specific 
values  a  and  y  for  o  and  y  that  minimize  J.  If  A  is  the  identity 
matrix,  then  the  minimization  is  in  the  least-squares  sense?  taking  A  as 
other  than  the  identity  matrix  allows  for  various  other  weightings  of 
the  difference  vector  Aph(o,y)  -  V.  A  standard  approach  for  such  mini¬ 
mization  problems  is  the  method  of  steepest  descent — a  method  of  itera¬ 
tively  refining  estimates  of  o  and  y  by  adjusting  their  values  dependent 
upon  the  behavior  of  the  gradient  of  the  cost  function.8 


8  A.  Sage  and  J.  Melsa,  System  Identification,  Academic  Press,  Inc., 
New  York  (1971). 
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The  method  revolves  critically  on  one's  ability  to  determine 
sequence  {k^}  of  constants  with  certain  convergence-insuring 
properties.  Specifically,  let  and  p^  denote  the  estimates  at  the  ith 
stage  of  the  descent,  and  let  =  J(a^,p^).  For  and  pi+1,  we  take 


(55) 


where  the  k^  must  be  such  that  Ji  1  Ji+1  at  every  stage  of  the  descent, 
and  equality  holds  only  in  the  case  ^  "  Ji+1  =  °*  The  explicit  con¬ 
struction  of  {  k^ }  for  a  given  problem  can  be  very  difficult  so  that  such 
sequences  are  often  established  by  empirical  investigation.  Such  inves¬ 
tigation  was  the  tack  chosen  for  the  problem  at  hand. 


Figure  17  shows  the  constant-cost  contours  of  J  for  A  =  the 
identity  matrix  and  for  the  true  values  o  =  0.2  m~*  and  p  =  0.01 
m_1sr-1.  Qualitatively  similar  pictures  are  obtained  for  other  true 
values  of  o  and  p.  Since  the  family  of  trajectories  orthogonal  to  the 
constant-cost  contours  gives  the  directions  of  the  vector  gradient  of  J, 
the  examination  of  such  pictures  can  suggest  appropriate  paths  for  the 
descent  and  point  out  potential  difficulties. 


The  steepest-descent  algorithm  for  k^  =  1/i  was  implemented 
on  a  computer  for  the  problem  with  true  values  a  =  0.2  m_1and  p  =  0.01 
m-1sr_1.  Trajectories  using  different  initial  step  sizes  (to  establish 
the  first  gradient  computations)  and  25  iterations  for  a  and  p  are  shown 
in  figure  18.  A  choice  for  the  matrix  A  that  weighted  the  values  of  the 
return  signal  more  heavily  near  its  peak  produced  similar  results. 


In  addition  to  steepest  descent,  an  intuitive  algorithm  for 
estimating  a  and  p  was  investigated.  The  algorithm  is  based  on  the 
observation  that  the  peak  of  the  return  signal  depends  more  strongly  on 
p,  while  the  decay  of  the  trailing  edge  depends  more  strongly  on  a. 
Accordingly,  the  following  algorithm  was  implemented  on  a  computer: 

a.  If  the  peak  signal  value  is  smaller  than  the  peak  value  of 
APh(o,p),  then  decrease  p;  otherwise,  increase  p. 


b.  If  the  trailing  edge  of  the  signal,  on  the  average,  decays 
more  rapidly  than  that  of  APh(o,p),  then  increase  a;  otherwise,  decrease 
a . 


The  sequence  { 1/i} was  used  again  to  scale  each  iteration.  The 
results  obtained  with  this  algorithm  jure  shown  for  several  typical 
computations  in  figure  19,  which  plots  trajectories  for  several 
different  step  sizes  and  for  the  same  true  values  of  a  and  p  used  in  the 
foregoing  examples. 
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Figure  17.  Constant-cost  contours  of  J(o,p)  = 

[Aph(o,n)  -  v] ' A [Aph( o , u)  -  V]  for  A  =  identity  matrix  and 
for  true  values  0  =  0.2  m-i  and  u  =  0.01  m-isr-i;  h(  a,  m)  is 
column  vector  with  components  tie~0cnR( cnA/2)  (n  =  0,  1, 

.  .  .,  99),  A  =  1  ns,  c  =  0.3  m/ns,  R  is  given  by  equation 
(29),  and  prime  denotes  transpose  operation.  Matrix  P  is 
defined  by  equation  (26)  with  P(t)  taken  as  in  figure  9. 
Fixed  devalues  are  shown  next  to  contours. 
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Figure  18.  Typical  trajectories  for  steepest 
descent  algorithm  with  ki  =  1/i.  True  values 
are  o  =  0.2  and  p  =  0.01  ln'^r"1,  and 
initial  guess  was  a  =  0.3  m-i  and  p  =  0.02 
m“Asr“A.  Several  initial  step  sizes  were  used 
and  in  each  case  25  iterations  of  descent  are 
shown. 


EXTINCTION  COEFFICIENT,  a  (m  -  >) 


Figure  19.  Typical  trajectories  for  intuitive 
algorithm  that  adjusts  a  and  y  at  each  step  in 
accordance  with  comparison  of  peak  values  and 
decay  times  of  APh( o,y)  and  true  signal.  True 
values  and  initial  guess  are  as  in  example  of 
figure  18.  Results  for  several  initial  step 
sizes  are  shown;  sequence  {l/i}  was  used  to 
scale  corresponding  iteration. 


5.2 


Estimation  of  a  and  tj  from  Measured  Return  Signals 


To  further  investigate  the  utility  of  parameter  estimation, 
the  steepest-descent  minimization  procedure  was  applied  to  10  samples  of 
measured  returns  from  cumulus  clouds.  These  sample  signals  were  taken 
from  the  HDL  data  collection  (sect.  1)  and  were  selected  to  correspond 
to  measurements  taken  with  the  GaAs  probe  fully  immersed  in  approxi¬ 
mately  uniform  cloud,  as  determined  by  another  instrument.  The  matrix  A 
defining  the  cost  function  J  of  equation  (54)  was  chosen  to  weight 
APh(a,y)  -  V  in  proportion  to  the  measured  signal  value.  This  choice 
makes  the  model  pulse  APh(a,u)  fit  the  measured  one  better  where  the 
signal  level  is  high,  thus  deemphasizing  the  effect  of  measurement 
noise.  Hie  sequence  {k^}  used  to  scale  the  descent  was  established 
empirically  and  has  the  form 

ki  =  c/(20  +  A)  '  (56) 

where  the  constant  is  determined  by  normalization  conditions. 

Figure  20  (a  to  h)  illustrates  the  results  obtained  for  eight 
of  the  sample  signals.  The  solid  curves  in  these  graphs  represent  APh 
for  the  best-fit  values  of  cr  and  u ,  and  the  dots  show  the  sampled  values 
of  the  measured  return.  The  results  for  the  remaining  two  sample 
signals  gave  negative  best-fit  values  for  a.  Closer  examination  of 
additional  data  characterizing  the  cloud  regions  that  produced  these 
signals  revealed  that  the  cloud  distributions  were  decidedly  nonuniform. 

To  determine  the  practicality  of  applying  the  foregoing 
approach  to  multilayer  models  needs  further  investigation.  While  there 
exist  many  computer  routines  for  determining  the  extrema  of  functions  of 
many  variables,  such  routines  are  typically  structured  around  a  specific 
problem  and  might  not  be  effective  for  other  problems.  An  investigation 
of  the  structure  of  the  M-layer  analog  of  equation  (54)  together  with  a 
review  of  available  minimization  routines  should  clarify  this  issue. 
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Figure  20.  Sampled  values  ot  measured  cumulus  cloud  return 
signal  (dots)  obtained  with  transmitter  pulse  of  figure  1 
and  range-response  curve  of  figure  3.  Solid  curves  give 
APh(a,p)  for  best-fit  values  of  o  and  p  determined  by 
steepest-descent  algorithm,  where  measurement  system  is 
assumed  to  have  been  fully  immersed  in  uniform  aerosol  of 
unknown  extinction  and  backscatter  coefficients. 
Appropriate  sampling  of  figure  1  transmitter  pulse  and 
equation  (29)  approximation  ot  figure  3  range  response  was 
used  for  algorithm  computations. 


Figure  20  (Cont'd).  Sampled  values  of  measured  cumulus 
cloud  return  signal  (dots)  obtained  with  transmitter  pulse 
of  figure  1  and  range-response  curve  of  figure  3.  Solid 
curves  give  APh (o,p)  for  best-fit  values  of  a  and  p 
determined  by  steepest-descent  algorithm,  where  measurement 
system  is  assumed  to  have  been  fully  immersed  in  uniform 
aerosol  of  unknown  extinction  and  backscatter  coefficients. 
Appropriate  sampling  of  figure  1  transmitter  pulse  and 
equation  (29)  approximation  of  figure  3  range  response  was 
used  for  algorithm  computations. 


Figure  20  (Cont'd).  Sampled  values  of  measured  cumulus 
cloud  return  signal  (dots)  obtained  with  transmitter  pulse 
of  figure  1  and  range-response  curve  of  figure  3.  Solid 
curves  give  APh(o,u)  for  best-fit  values  of  o  and  p 
determined  by  steepest-descent  algorithm,  where  measurement 
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Figure  20  (Cont'd).  Sampled  values  ol  measured  cumulus 
cloud  return  signal  (dots)  obtained  with  transmitter  pulse 
of  figure  1  and  range-response  curve  of  figure  3.  Solid 
curves  give  APh(o,p)  for  best-fit  values  of  o  and  p 
determined  by  steepest-descent  algorithm,  where  measurement 
system  is  assumed  to  have  been  fully  immersed  in  uniform 
aerosol  of  unknown  extinction  and  backscatter  coefficients. 
Appropriate  sampling  of  figure  1  transmitter  pulse  and 
equation  (29)  approximation  of  figure  3  range  response  was 
used  for  algorithm  computations. 
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Figure  20  (Cont’d).  Sampled  values  of  measured  cumulus 
cloud  return  signal  (dots)  obtained  with  transmitter  pulse 
of  figure  l  and  range-response  curve  of  figure  3.  Solid 
curves  give  &Ph(o,u)  for  best-fit  values  of  o  and  m 
determined  by  steepest-descent  algorithm,  where  measurement 
system  is  assumed  to  have  been  fully  immersed  in  uniform 
aerosol  of  unknown  extinction  and  backscatter  coefficients. 
Appropriate  sampling  of  figure  1  transmitter  pulse  and 
equation  (29)  approximation  of  figure  3  range  response  was 
used  for  algorithm  computations. 
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Figure  20  (Cont'd).  Sampled  values  of  measured  cumulus 
cloud  return  signal  (dots)  obtained  with  transmitter  pulse 
of  figure  1  and  range-response  curve  of  figure  3.  Solid 
curves  give  APh(a,p)  for  best-fit  values  of  o  and  p 
determined  by  steepest-descent  algorithm,  where  measurement 
system  is  assumed  to  have  been  fully  immersed  in  uniform 
aerosol  of  unknown  extinction  and  backscatter  coefficients. 
Appropriate  sampling  of  figure  1  transmitter  pulse  and 
equation  (29)  approximation  of  figure  3  range  response  was 
used  for  algorithm  computations. 
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6.  SUMMARY  AND  DISCUSSION 


This  report  discusses  how  a  growing  HDL  data  bank,  consisting  ot 
measured  backscattered  laser  pulses  from  aerosols  such  as  weather  clouds 
and  smoke,  must  be  treated  before  such  data  can  be  directly  applied  to 
evaluate  the  aerosol  vulnerability  of  various  AOF  systems  and  methods 
for  aerosol  discrimination.  The  measurement  system,  a  short-pulse  GaAs 
laser  probe,  essentially  convolves  the  features  of  the  aerosol  that  are 
sought,  namely,  the  aerosol  signature,  with  the  shape  of  the  probe’s 
transmitter  pulse,  after  scaling  the  signature  by  the  system's  range- 
response  characteristic.  Thus  deconvolution  and  rescaling  are  necessary 
to  obtain  the  desired  aerosol  signatures. 

Rescaling  of  the  data  using  the  known  range-response  characteristic 
of  the  probe  is  straightforward.  Deconvolution  with  respect  to  the 
transmitter  pulse,  while  possible  in  principle,  is  difficult  to  accom¬ 
plish  in  practice,  unless  one  is  willing  to  accept  the  spatial  resolu¬ 
tion  implied  by  the  width  of  the  probing  pulses  (*  2  m).  In  the  latter 
case,  deconvolution  reduces  to  a  mere  additional  rescaling. 

Two-meter  spatial  resolution  is  acceptable  for  most  purposes,  but  is 
a  serious  drawback  in  attempting  to  gauge  the  buildup  rates  of  aerosol 
density  near  the  boundaries  of  an  aerosol  distribution.  The  deleterious 
effects  of  sharply  rising  aerosol  edges  on  AOF  system  performance  are 
likely  the  most  difficult  ones  to  cope  with  in  designing  an  aerosol 
resistant  system.  It  is  therefore  of  great  interest  to  obtain  high- 
resolution  measurements  of  aerosol  signatures  in  regions  where  they  may 
be  changing  rapidly. 

In  accepting  the  resolution  implied  by  the  width  of  the  transmitter 
pulse,  one  assumes  as  an  approximation  that  the  transmitter  pulse  can  be 
considered  as  a  spatial  5-function  for  the  purpose  of  unraveling  its 
convolution  with  the  range-response  scaled  aerosol  signature.  This 
approximation  results  in  the  interpretation  of  the  backscattered  pulse 
as  being  proportional  to  the  scaled  aerosol  signature,  and  in  a  roughly 
2-m-resolution  measurement.  It  is  not  necessary  to  make  the  foregoing 
approximate  assumption,  however.  If  the  actual  shape  of  the  probing 
pulse  is  taken  into  account,  the  resulting  convolution  relation  between 
the  backscattered  pulse  and  the  aerosol  signature  not  only  is  definite, 
but  also  is  solvable  for  the  signature,  at  least  in  principle.  Conse¬ 
quently,  increased  resolution  is  theoretically  possible  through  a  more 
sophisticated  interpretation  of  the  data.  The  bulk  of  this  report 
concerns  the  investigation  and  the  development  of  this  more  sophisti¬ 
cated  view.  The  findings  and  the  results  are  summarized  below. 

It  has  been  demonstrated  that  the  mathematical  process  required  to 
deconvolve  backscattered  pulses  can  be  economically  implemented  on  a 
computer.  An  algorithm  that  uses  discrete  Z-transform  methods  and  is 
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ideally  suited  to  the  digital  form  of  the  backscatter  dati  has  been 
coded  on  computer  and  shown  to  be  effective.  Moreover,  the  running  time 
of  the  program  is  sufficiently  short  to  make  the  deconvolution  of  a 
large  data  bank  feasible. 

The  essential  difficulty  in  determining  aerosol  signatures  by  decon¬ 
volution  is  the  effect  of  the  noise  present  in  the  backscattered  pulse 
measurements.  A  complete  analysis  of  the  effects  of  noise  has  been 
carried  out  and  shows  that  the  magnitude  of  noise  errors  in  the  aerosol 
signature  depends  sensitively  on  the  noise  spectrum  in  relation  to  the 
Fourier  transform  of  the  transmitter  pulse.  For  example,  if  the  latter 
vanishes  at  certain  frequencies,  and  if  the  backscattered  pulse's  noise 
spectrum  does  not  vanish  sufficiently  ’'epj.dly  as  these  frequencies  are 
approached,  then  the  mean-squared  noise  level  associated  with  the 
aerosol  signature  becomes  infinite;  that  is,  the  SNR  of  the  decon¬ 
volution-determined  aerosol  signature  is  zero.  This  catastrophic  result 
can  be  avoided  by  appropriate  filtering  of  the  backscatter  data  prior  to 
performing  Jthe  deconvolution  and  does  not  occur  at  all  if  the  Fourier 
transform,  P(f),  of  the  transmitter  pulse  does  not  vanish  within  the 
bandpass  of  the  measurement  system  receivers.  Moreover,  P(f)  does  not 
appear  to  vanish  within  the  relevant  bandpass,  which  is  approximately 
200  MHz,  since  a  numerical  Fourier  transform  of  a  typical  measured 
transmitter  pulse  (fig.  2)  showed  no  zeros  for  frequencies  up  to  250 
MHz.  However,  backscatter  data  have  been  obtained  with  several  distinct 
transmitter  pulse  shapes  so  that  further  analysis  is  required  on  this 
point. 

The  noise-effect  analysis  has  been  applied  to  the  Fourier  spectrum 
exhibited  by  the  measured  transmitter  pulse  of  figure  2.  The  relation¬ 
ship  between  the  SNR  (peak  signal  divided  by  rms  noise)  of  a  back- 
scattered  pulse  and  that  of  the  corresponding  deconvolution-determined 
aerosol  signature  has  thereby  been  determined  (fig.  16).  For  the 
approximately  200-MHz  receiver  bandwidth,  the  results  indicate  that  the 
SNR  of  the  aerosol  signature  is  0. 15  times  the  SNR  of  the  backscattered 
pulse.  To  obtain  reasonably  good  signature  SNR's  in  the  face  of  such  a 
reduction  factor  requires  a  very  good  SNR  for  the  backscatter  measure¬ 
ment.  Since  the  typical  SNR  for  the  backscatter  measurements  obtained 
with  the  transmitter  pulse  being  considered  is  only  about  15  to  1,  we 
must  conclude  that  the  deconvolution  method  will  fail  to  give  acceptable 
results  for  these  data . 

Smoothing  of  the  data  has  been  considered  as  a  way  of  improving  the 
signature  SNR;  however,  smoothing  is  tantamount  to  reducing  the  measure¬ 
ment  bandwidth.  The  receiver  bandwidth,  which  is  actually  somewhat 
greater  than  200  MHz,  results  in  a  smearing  of  the  received  optical 
pulses  by  a  spatial  average  over  about  .a  0.5-m  interval.  We  term  this 
averaging  interval  the  bandwidth  resolution.  Since  2-m  resolution  of 
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the  aerosol  signatures  is  already  available  without  detailed  decon¬ 
volution,  there  is  no  advantage  in  smoothing  to  a  bandwidth  resolution 
approaching  2  m.  The  effect  of  smoothing  for  the  case  in  figure  16  can 
be  seen  by  considering  what  would  occur  if  the  bandwidth  were  reduced  to 
150  MHz.  Then  the  bandwidth  resolution  would  worsen  to  about  1  m,  while 
the  signature  SNR  would  be  improved  by  only  about  30  percent  (fig.  16). 

The  most  recent  series  of  backscatter  measurements  has  been  made  by 
using  a  shorter  (5-ns  FWHM) ,  more  symmetrical,  and  considerably  higher 
peak-power  transmitter  pulse  than  that  shown  in  figure  2.  Typical  SNR's 
for  the  backscattered  pulses,  while  not  yet  determined  accurately  from 
the  data,  are  expected  to  be  5  to  10  times  greater  than  for  previous 
measurements  and  therefore  should  make  the  deconvolution  method  more 
attractive.  A  noise  analysis  using  the  new  transmitter  pulse  must  still 
be  made;  but  if  the  results  are  not  too  different  from  those  for  the 
pulse  of  figure  2,  then  signature  SNR's  approaching  20  to  1  may  be 
obtainable . 

Better  signature  extraction  for  the  older  data  may  be  possible 
through  an  alternative  approach.  A  preliminary  investigation  of  using  a 
priori  knowledge  about  the  signature — mainly  its  definition  in  terms  of 
the  physical  parameters  of  the  aerosol — and  parameter  estimation  tech¬ 
niques  has  shown  that  exellent  results  can  be  obtained  in  simple  cases 
(where  the  aerosol  is  known  to  be  uniformly  distributed  and  to 
completely  engulf  the  measurement  system) .  The  method  applies  to  all  of 
the  backscatter  data  presently  available  and,  in  addition,  seems  capable 
of  providing  more  detailed  information  about  the  aerosol  than  is  con¬ 
tained  in  its  signature.  The  computational  aspects  of  the  method, 
however ,  are  considerably  more  complex  than  those  of  direct  decon¬ 
volution  and  may  be  difficult  to  implement  for  highly  nonuniform 
aerosols. 
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